Draft version November 29, 2010 

Preprint typeset using IAT^X style cmulateapj v. 11/10/09 



CHARACTERIZING THE TIME VARIABILITY IN MAGNETIZED NEUTRINO-COOLED ACCRETION 
DISKS: SIGNATURES OF THE GAMMA-RAY BURST CENTRAL ENGINE 

Augusto Carballido and William H. Lee 

Instituto de Astronomi'a, Univcrsidad Nacional Autonoma de Mexico 
Apdo. Postal 70-543, Mexico D.F. 04510, MEXICO 
Draft version November 29, 2010 

ABSTRACT 

The central engine of Gamma Ray Bursts is hidden from direct probing with photons mainly due 
to the high densities involved. Inferences on their properties are thus made from their cosmological 
setting, energetics, low-energy counterparts and variability. If GRBs are powered by hypercritical 
accretion onto compact objects, on small spatial scales the flow will exhibit fluctuations, which could 
in principle be reflected in the power output of the central engine and ultimately in the high energy 
prompt emission. Here we address this issue by characterizing the variability in neutrino cooled 
accretion flows through local shearing box simulations with magnetic fields, and then convolving 
them on a global scale with large scale dynamical simulations of accretion disks. The resulting 
signature is characteristic, and sensitive to the details of the cooling mechanism, providing in principle 
a discriminant for GRB central engine properties. 

Subject headings: accretion, accretion disks — gamma rays bursts: general — hydrodynamics - 
magnetic fields 
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1. INTRODUCTION 

Ray Bursts (GRBs) are thought 
driven by mass accretio n onto 
objects dWooslev fc Blooml 
iLee fe Ramirez-Ruia 



obje 
120(171 



2006 



2007 



Gehrels. Ramirez-Ruiz fe Foxll2009l ). The enabling cool- 



ing mechanism allowing accretion is neutrino emission, 
raising the usual Eddi ngton rate for photons by sixteen 
orders of magnitude (IPopham. Wooslev fe Frverl 119991: 
Naravan. Piran fc Kuman 120011: iKohri fc Mincshiac 
2002t iBeloborodovl 120031: iDi Matteo. Perna fc Naravan! 
20021 iSetiawan. Ruffert fc Jankal 120041: 

Lee. Ramirez-Ruiz fc Page] 120051 : iChen fc Beloborodovl 
20071) . This applies both to long and short GRBs, 
erhaps due to t he collapse of massive r o tating stars 



pernaps d ue to tr ie collapse ol massive r< 
(|Wooslevl 119931 : iMacFadven fc Wooslevi 



.19991) and 

compact binary mergers (lEichler et alJll989t " p*aczvnskil 
119911 : iNarayan. Paczynski fc Piranl I1992T) . respectively, 
or magnetized neutron stars (|Usovill992l ). The prompt 
gamma-ray emission originates at ~ 10 14 — 10 16 cm from 
this source, possibly from internal sh ocks in the relativis- 
tic ou tflow generated by the engine (IZhang fc MeszarosI 
2004). The engine itself is hidden from view due to the 
high opacities, and can be potentially probed directly 
only through gravitational waves and neutrinos. 

The observed time series in GRB pro mpt emission 
show diversity between events (see, e.g., iNorris et al.l 
[1996|): some have a single peak, others multiple emission 
episodes, with correlations between the fluence of the 
active pe riod and the length of the quiesc ent interval pre- 
ceding it (Ra mirez-Ruiz fc Merlonil2001l) . On top of this, 
rapid (ms) fluctuations are routinely observed. Fourier 
analysis of the high-energy light curves of the prompt 
emission i n GRBs in the source frame reveals po w er-law 
spectra ([B elobor odov. Stern fc Svenssonl 119981 120001: 
iRvde et al.ll2003D . with index ~ —5/3 and a break at 
~ 1 — 2 Hz. Variability is likely due to a combination of 



several effects, allowing in principle an additional way to 
probe the central engine indirectly. Some are probably 
intrinsic to the progenitor: the distribution of angular 
momentum with radius inside the star in the case of a 
collapsar may lead to distinct episodes of energy release 
(iKuma r. Nara van fc Johnson 2008; Perna fc MacFadvenl 
120101 : lLopez-Camara. Lee fc Ramirez-Ruizl I2010D : 
the fall back at late times of material stripped 
from a tidally disrupted neutron star is capable 
of powe ring secondary accretion episodes dRosswogl 



120071 ILee. Ramirez-Ruiz fc Lopez-Camaral 120091) : 
hydrodynamical or magnetic instabilities in the ac- 
cretion disk may result in i ntermittent accretion 
(iPerna. Armitage fc Zhangl 120061: iPrqga fc Zhangl 120061 : 
iTavlor. Miller fc Podsiadlowskil boiOt ) Others can 
come from the relativistic outflow: the interaction 
of a jet with high Lorcntz factor with the stellar 
envelope before breakout can lead to irreg u laritie s 
and shocking (|Morsonv. Lazzati fc Begelmanl [2010); 
the outcome of internal shocks between shells in 
the flow depends on t he variation in mass and en- 
ergy upon ejection dPanaitescu. Spad a & M eszaroi 
1999; iRamirez-Ruiz. Merloni fc Reesl 120011 : 



Bosniak. Daigne fc Dubusl 120091: iMendoza et al.l 1200 
An additional factor, upon which we focus here, is 
related to the variability present in the accretion disk as 
a result of turbulent motions. The dissipation is related 
to the local hydrodynamical variables, and as these vary 
in time, so will the energy output. 

The magnet o rotati onal instability (MRI) 
(Bal bus fc Hawlevl fl998h is a possible mechanism 
that will allow for the transport of angular momen- 
tum in accretion disks with differential rotation. Its 
behavior under the physical conditions in neutrino 
cooled disks, which are similar to those occurring in 
supernovae, but different from the usual ones present 
in X-ray binari es and AGN, has come under scrutiny 
more recently ([Thompson. Quataert fc Burrows! 120051 : 
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Masada et al.l 120071: iRossi. Armitage & Meiioul 120081: 
Obergaulinger et al.l I2009D . One such difference lies in 
the sensitivity of the cooling rate to temperature and 
is at the heart of this work. Whereas, for example, the 
photon bremsstrahlung emissivity scales as q oc T 1 / 2 in 
the optically thin limit, for neutrinos q oc T* 3 , where 
P ~ 6 — 9 depending on the cooling process. Further, 
while photon-cooled disks are typically optically thick, 
t 7 3> 1, for a wide range of relevant parameters their 
neutrino-cooled counterparts are optically thin, t u < 1, 
tightly coupling the local conditions to the emitted 
luminosity. 

In this Letter, we characterize the local variability in 
neutrino cooled disks through shearing box MHD sim- 
ulations (§ 0), and then use the results of global disk 
simulations to scale the results for the central engine as 
a whole ( § EJ) - Prospects for placing constraints on GRBs 
are discussed in § 111 

2. THE SMALL SCALE: NUMERICS AND PHYSICS IN THE 
SHEARING BOX CALCULATIONS 

The local flow in the disk is modeled by the shear ing 
box approximation (Ha wlev. Gammie fc B albus 1995|): a 
rectangular Cartesian coordinate system represents a lo- 
cal neighborhood inside the disk, at an arbitrary orbital 
radius Rq with dimensions which are much smaller than 
Rq. The x, y and z axes represent the radial, azimuthal 
and vertical directions in the disk, respectively. The ra- 
dial component of the central object's gravity is included, 
and differential rotation is replaced by a Keplerian shear 
flow along the y direction. Periodic boundary conditions 
are set in y and z, while in x these are "shearing peri- 
odic" : upon crossing a radial boundary, a fluid element 
is displaced along y by an amount gi yen by the shear 
value. The Eulerian ZEUS-3D code [Stone k Normanl 
(1992a. b)] is used to solve the equations of ideal MHD. 
These are: 

^ + V-(pv)=0, (1) 



T 6 . A fraction £ of the initial internal energy of the fluid 
is removed over 100 orbital periods at Ro, the duration 
of the simulation. We used £ = 0.1, thus the cooling time 
is much longer than the dynamical time. 

Initially the density, internal energy and pressure in 
the box are uniform, and the magnetic field is vertical 
with non-zero net magnetic flux across the box. The 
ratio of gas to magnetic pressure is P ga s/-Pmag = 400. 
The velocity field is randomly perturbed with amplitude 
Wpcrt = 10~ 3 c s , where c s is the sound speed. We have 
performed simulations at varying resolution to test for 
convergence (the highest having 128x256x128 zones in 
x,y,z, respectively). From the power spectrum of the 
velocity field Vi(ki), we have verified that our results are 
converged and unaffected by resolution changes at the 
level reported. Turbulent motions induced by the MRI 
develop in the simulation on a time scale £mri cx l/f^Oi 
leading to turbulent transport of angular momentum. 
The fundamental output of the simulations is the lumi- 
nosity as a function of time of the shearing box, which 
we denote as s(t). 
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Fig. 1. — Left: Energy output (arbitrary units) as a function of 
time, s(t), for shearing box simulations with cooling mimicking pair 
capture (blue, /3 = 6) and pair annihilation (red, /3 = 9). The time 
is given in units of the orbital period at Ro, the fiducial shearing 
box radial position. Right: Fourier transforms s(u) of s(t). The 
reference power law has index -2. 



<9v _ 1 / B 2 \ (B-V)B _ 3^ 2 

(2) 

— = Vx(vxB), (3) 

g = -PVv-C^, (4) 

and an ideal gas equation of state with P — (7— l)e where 
v, p and P are the gas velocity, density and pressure, re- 
spectively, B is the magnetic field, fio = (0, 0, SIq) is the 
disk angular frequency, e is the internal energy density, 
and x is a unitary radial vector. 

The second term on the right-hand side of Eq.(j4]) is 
added to model neutrino cooling processes, with the 
nature of the particular emission mechanism fixing the 
value of p. In GRB central engi nes, the dominan t pro- 
cesses are e ± pair annihilation (Ito h et al.l 11996). with 
(7 ann oc T 9 , and capture onto free neu trons and pro- 
tons (|Langanke fc Martinez-Pinedoll200iT ). with g cap oc 



Luminosity variations for s(t) are shown in Figure Q] 
'for j3 = 6,9. Rapid variability is present, reflecting the 
turbulent motions in the box. These changes can be char- 
acterized by their Fourier transform, s(v), where v is the 
signal frequency (Figured]). Roughly, the spectrum is a 
power law, s(v) oc v~ 2 , over a wide range of frequencies 
in the case of /? = 6, with more complicated behavior su- 
perimposed in the case of (5 = 9. It is clear that the vari- 
ations in energy release have a characteristic behavior, 
with the details varying according to the specific cooling 
mechanism. 

3. THE LARGE SCALE: EXTENSION TO GLOBAL DISK 
VARIABILITY 

To find the global variability we require the superposi- 
tion of a very large number of smaller regions represented 
by the shearing box calculations, and approximations are 
necessary to proceed further. We assume equatorial and 
azimuthal symmetry, and also take the disk to be in a sta- 
tionary state. Obtaining global quantities thus requires 
integrating over the polar angle S [0, 27r] and radius 
r £ [ri n ,oo]. We are still faced with the superposition 
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of a potentially large number of zones, each uncorrelated 
with the other, and which could destroy any coherent 
signature of the processes occuring on small scales. A 
natural quantity in the disk, however, enforces causality 
for a given frequency: the local sound speed, c s . 

Consider a disk annulus at radius R, with sound 
speed c s . The correlation length over which signals with 
frequency v propagate is l v = c s jv. The azimuthal 
velocity is = and from vertical hydrostatic bal- 
ance, the pressure scale height H satisfies H/R ~ Cg/v^. 
Thus c s = 2irRf(H/R), with / being the orbital 
frequency at radius R. In neutrino cooled accretion 
disks, A = H/R est (ILee. Ramirez- Ruiz fc Pagdl2005t 
iLee fc Ramirez-Ruidl2007l ). so that L = 2irRAflv. The 
number of distinct azimuthal zones within the ring for 
frequency v is N$ lV = 2ttR/I l , — v/(Af). Dynamical 
simulations of such disks show that A » 1/10. If 
Nq.v > 1) the signal at frequency v is thus diminished 
with respect to the idealized case of entirely in-phase 

1/2 

contributions by a factor NJ V , due to the random nature 
of the phase de-correlation. If, on the contrary, N^, v < 1, 
the signal at frequency v is added constructively. This 
acts as a broad low-pass filter by slightly suppressing the 
signal at high frequencies, when v ^> Af. Carrying out 
this filtering procedure also requires scaling the fiducial 
model to different radii, which we now address. 

The emissivity is not a function of azimuthal position, 
but it is one of radius. Given the extreme sensitivity of 
the cooling rate to temperature, only the innermost re- 
gions of the disk contribute to the luminosity, and thus 
to its variability properties. The argument given above 
for the decoupling of azimuthal zones in the disk for a 
particular frequency is also applicable to the radial coor- 
dinate. Advection may modify this slightly, but for our 
present purpose it is valid enough to say that as the flow 
is split up into zones in azimuth, it is divided into 
rings of characteristic size l v . The number of radial zones 
in the disk at frequency v is Nr v = r ou t/l v , where r out is 
the outer radius of the disk. Using the previous results, 
Nrv — r ou tV / '{2-KRAf). A second low- pass filter needs 
to be considered in analogy with the azimuthal one, and 

1 /2 

the signal is now suppressed by a factor N^ v if Nr v > 1. 

In order to carry out this procedure, we have used 
the results of global , 2D simulations of neutrino- 
cooled accretion disks (ILee. Ramirez-Ruiz fc Pagdl2005t 
ILee fc Ramirez- Ruiz! 1200^ ) orbiting a stellar mass black 
hole, evolved in t he p seudo-relativistic potential of 
iPaczvhski fc W iita (198(3). For consistency purposes, we 
have re-computed their evolution for this work for two 
test cases, with Mbh = 4M Q and q = M^s^/Mbh = 
7.5 x 10~ 3 . In the first model, T6, we have eliminated 
all sources of cooling except from pair captures, where 
q ocT 6 , while in the second, T9, we have only kept that 
due to pair annihilation, with q oc T 9 . Each model can 
then be more faithfully compared with the two fiducial 
runs for the shearing box with j3 = 6, 9 respectively, us- 
ing a near-stationary state where the accretion rate is 
M = 0.05 Mq s _1 . 

We first compute the equatorial radial profile of emis- 
sivity in simulations T6 and T9, q(R), shown in Fig- 
ure [2] At small radii the dominant contribution comes 
from pair capture, while annihilation dominates for R > 



1e+32 




Fig. 2. — Equatorial [z = 0) neutrino emissivities, q(R/Rsch)i 
Rsch = 2GMbh/c 2 , for two dimensional, azimuthally symmetric 
global dynamical simulations of accretion disks around a black hole 
with M = 4Mq. The blue (red) line corresponds to simulation T6 
(T9), where only pair captures (annihilation) were used to compute 
the cooling rate. The solid black line is the result of a simulation 
including both processes. 

R* = 30i?s c h. The combination is a strongly decreas- 
ing function of R, so the radial convolution need not be 
extended to infinity, and pair capture dominates the to- 
tal emission, as well as its general behavior. The inner 
radius of the disk is r- m = r ms = 3i?sch, where the ef- 
fects of General Relativity truncate the flow. We have 
used r out = 85i?sch, and checked that the final result is 
insensitive to this choice, as long as r out ^> R*. 

To perform the radial integral, the time series s(t) and 
s{v) need to be scaled to the disk at different radii R. The 
characteristic break frequency in the fiducial spectrum 
for the volume-averaged total (Maxwell and Reynolds) 
stress, measured as an equivalent a viscosity parameter, 
is roughly at the orbital frequency, /. This is related 
to the nature of the MRI, where, £mri — I/O,. The or- 
bital angular frequency is known as a function of radius, 
SI = (GM B n/R) 1/2 (R/{R - Rsch)), so we can convert 
from frequency to radius, s(v) to s(R) and integrate over 
different annuli (the expression for Q is appropriate for 
the pseudo-relativistic potential used). 

The power associated to the variability in the neutrino 
luminosity can now be computed as: 

P{v) = f q(R)3>(N (t>u )p(N Rv )s(R)2TrRdR, (5) 

J frill 

where q(R) corresponds to model T6 or T9, and ^(N^^) 
and p(Nn u ) account for the filtering at high frequen- 
cies. Figure [3] shows the resulting spectra, where sev- 
eral features clearly stand out. At low frequencies, 
v < v o ~0.1 Hz, the spectrum is the same in both cases, 
P{v) oc v~ . For v > vq the results diverge: model 
T6 approximately maintains the power law decay, while 
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Fig. 3. — Convolved power spectra over the entire disk. Eq. [5] 
for models T6 (red) and T9 (blue). The power is given in arbitrary 
units as a function of frequency (in Hz) for a black hole with 4Mq . 
The black line shows the result of a test calculation with local 
cooling in the shearing box simulation set to q oc T, i.e., /3 = 1. 
The error bars for = 6, 9 include the scatter from the Fourier 
transform in the plotted frequency bins and sampling errors in the 
radial discretization of the convolution. 

model T9 shows an excess, and in general more complex 
behavior, with broad features at ~ 2 and 40 Hz overlaid 
on a decay with lower index, P(v) oc v~ lr ' ' . At v = V\ ~ 
100 Hz, a break to a more rapidly decaying power law ter- 
minates both spectra. The errors are greater for model 
T9, due to greater scatter in the Fourier transform. The 
thin black line in Figure [3] is from a shearing box test 
calculation with q oc T, i.e., with /3 = 1 (the smaller er- 
ror bars are not plotted for clarity). As the sensitivity 
of the cooling term to the temperature rises, the trend is 
for an excess in power at higher frequencies, and a more 
moderate background power law decay. 

4. SUMMARY, DISCUSSION AND PROSPECTS FOR 
OBSERVABILITY IN GRBS 

Through shearing box MHD simulations we have char- 
acterized the time variability of the energy release at 
small scales in a neutrino-cooled accretion disk around 
a black hole, where energy losses primordially come from 
e^- pair capture onto free nucleons and protons and pair 
annihilation. With the use of cooling profiles from large 
scale, two-dimensional simulations of full disks, we have 
convolved this local variability to obtain a global signa- 
ture of time variations in the power output, through the 
power density spectrum of the neutrino luminosity. Since 
accretion is enabled by the cooling through neutrinos, 
we take this as an indicator of central engine variability 
which will be reflected in the relativistic outflow eventu- 
ally giving rise to a GRB. 

The power spectrum exhibits characteristic features re- 
lated to the general nature of the neutrino cooling in the 
optically thin regime, and particularly to its temperature 



dependence. A background power law decay, with index 
~ 1.7 — 2 extends approximately from 0.1-100 Hz. For 
emissivities q oc T 13 , additional power appears at high 
frequencies as j3 rises, with the particular values of the 
slopes and cuttoffs scaling with the cooling mechanism 
(Figj3|). This result thus provides in principle a discrimi- 
nant for GRB central engines powered by neutrino-cooled 
accretion flows, and illustrates how any local mechanism 
can be used in the same way to test its viability. 

The reason for the difference when modifying the cool- 
ing mechanism is fundamentally related to the local en- 
ergy balance. For q oc T 9 (model T9) and a given inter- 
nal energy supply, e, the local cooling time t coo \ — e/q is 
shorter under a temperature perturbation AT than when 
q oc T e (model T6) . The power is thus higher at such fre- 
quencies, leading to the observed displacement in Fig. |3l 
The argument holds when comparing model T6 with the 
test run at /3 = 1 , and gives a way to characterize the ac- 
cretion flow and discriminate between competing mech- 
anisms, if the variations in the neutrino luminosity are 
directly reflected in the accretion power output which 
drives a relativistic flow. 

Two possible scenarios in which the difference between 
the cooling regimes studied here may be of relevance are 
worthy of note. First, the mass of the black hole intro- 
duces a scaling into the problem when combined with 
neutrino cooling. If Mbh is too large, the density and 
temperature in the accretion flow can be too low for neu- 
trino cooling to operate. As the BH mass is reduced, pair 
annihilation, with q oc T 9 first becomes effective as an 
energy sink, followed by pair captures, with q oc T 6 . A 
signature of the BH mass is thus in principle available in 
the variability of the flow. The second case is related to 
the time evolution of the flow, assuming a certain amount 
of mass Mdi s k is initially available for accretion, and no 
further feeding of the central engine takes place. Global 
disk simulations show that the density and temperature 
drop as the disk drains into the BH on the viscous time 
scale. Given enough time, the whole disk will lie on the 
branch cooled by pair annihilation where a particles have 
formed, below log[p(g cm -3 )] ~ 6.5 and log[T(°K)] ~ 10 
for the adopted black hole mass (see also Figure 2 in 
iLee. Ramirez-Ruiz k. Lopez-Camaral (120090 1. Thus the 
variability may initially behave as in case T6, and end 
as in case T9 (the luminosity will have decreased sub- 
stantially by then, along with the accretion rate). The 
associated GRB, if one occurs, need not necessarily be 
powered by neutrinos in order for these effects to be ap- 
parent. The fact that neutrinos are responsible for the 
cooling process allowing accretion makes them relevant 
in this context. 

A number of limitations apply to this study, and can be 
matters for further study. First, the scaling we have used 
to infer the variability at different radii is strictly valid if 
the flow is adiabatic. Explicit cooling thus violates this 
assumption. However, given the choice of £ and i coo i the 
associated time scales are such that i coo i 3> idyn, making 
this a reasonable approximation. Second, we have used 
the neutrino luminosity as a proxy for the manifestation 
of central engine activity, which could have a neutrino 
component, but need not be restricted to it. For ex- 
ample, magnetic fields may power relativistic outflows 
leading to the observed high-energy emission. In this 
sense, an alternative characterization based on the mass 
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accretion rate M through the disk may be of use as well, 
which can then be associated to the energy output of the 
relativistic outflow as L TC \ oc Mc 2 . Third, this is only the 
first filter any variability originating in a disk needs to go 
through before emerging as high-energy photons. A vari- 
able mass accretion rate and energy conversion efficiency, 
propagation through the stellar envelope (for a massive 
star pr ogenitor), externa l shocks, and general relativistic 
effects (|Birkl et al.ll2007tl . among others, can each poten- 
tially leave their own fingerprint on the power spectrum. 
The background upon which they will do so, however, 
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